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Abstract 

We present exact methods for the numerical integration of the Wannier-Stark system in a many-body sce¬ 
nario including two Bloch bands. Our ab initio approaches allow for the treatment of a few-body problem 
with bosonic statistics and strong interparticle interaction. The numerical implementation is based on the 
Lanczos algorithm for the diagonalization of large, but sparse symmetric Floquet matrices. We analyze the 
scheme efficiency in terms of the computational time, which is shown to scale polynomially with the size 
of the system. The numerically computed eigensystem is applied to the analysis of the Floquet Hamiltonian 
describing our problem. We show that this allows, for instance, for the efficient detection and characteri¬ 
zation of avoided crossings and their - statistical analysis. We finally compare the efficiency of our Lanczos 
diagonalization for computing the temporal evolution of our many-body system with an explicit fourth or¬ 
der Runge-Kutta integration. Both implementations heavily exploit efficient matrix-vector multiplication 
schemes. Our results should permit an extrapolation of the applicability of exact methods to increasing 
sizes of generic many-body quantum problems with bosonic statistics. 

Keywords: Floquet problem, Bose-Hubbard model, Wannier-Stark system, Lanczos algorithm. Quantum 
chaos, Ultracold atoms. Optical lattices. 
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1. Introduction 

1.1. The Wannier-Stark problem 

The problem of a particle in a periodic potential and subjected to a constant field has a long history in 
quantum physics. Today known as Wannier-Stark problem HHH, it is the basic model for conductance. 
An appropriate matching between the field and the periodic lattice can induce resonance phenomena, and 
hence be used for the control of transport. This idea goes back to Esaki and Tsu who suggested to implement 
diodes and other transport devices with the techniques of semiconductor physics |3] |4j]. While of interest for 
transport, it was actually shown quite early that the constant field, also called tilt, induces regular periodic 
oscillations instead of an acceleration along macroscopic sizes. This phenomenon of Bloch oscillations and 
the fact that the Wannier-Stark problem is strictly speaking an open system has evoked much interest over 
decades fT|. 

Modern realizations of the Wannier-Stark system use either light |5| or ultracold atoms for which the 
potentials can be controlled very well f6j|7]], avoiding intrinsic effects of unwanted disorder and uncontrolled 
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particle-particle interactions in solids. Implementations with cold atoms have turned the attention from 
fermions also to bosonic transport particles J6][3, from single-particle models (2) over mean-field effects 
(U |9j |T0l im towards many-body correlation effects Ifl2l [T3l [T4l IT5l IT6l [171 [T8l IT9l l20ll . Needless to say 
that the many-body scenario is much more complicated. One typically has to resort then to effective models 
mm, which themselves build on a good intuition of the underlying physical phenomena that take place 
during the temporal evolution. The coupling between energy bands (Bloch bands in the periodic lattice) 
induced by the tilt makes effective single-body descriptions already non-trivial |ffi 2,22. [231 . not to speak 
about the growing size of the system when including many particles in various dimensions. Here the Hilbert 
space grows exponentially, both with the size of the system and the number of particles. 

1.2. Our approach 

In view of the complications due to many coupled energy bands, spatial dimensions and many-body 
effects, we restrict here to a one-dimensional realization including fully two energy bands. Our system 
is thus closed, but it represents a good approximation to a possible realization with ultracold atoms in a 
double-periodic potential producing a miniband structure, see Fig. [T] below and the appendix in I! 171 . Here 
we present a detailed study of our computational methods which allow for strong interactions between the 
particles. We study finite systems, with a finite number of particles, bosons to be specific, in a finite number 
of sufficiently deep lattice sites. All these assumptions make it possible to use a lattice model, whereby 
the standard Bose Hubbard model Ell [24 1 is much extended here by including two coupled energy bands, 
interband interactions terms, and field-induced couplings. The exponential growth of the size of the problem 
with, e.g. the particle number, is especially related to the choice of bosons, which can, in principle, all sit 
on the same lattice site. This makes our problem non-integrable and computationally hard even in one 
spatial dimension lfl5l 47 1. Applying periodic boundary conditions, the constant field is mapped onto a 
periodically time-dependent quasi-momentum, making our system effectively driven by a periodic force. 
This corresponds to an acceleration of the lattice as done also experimentally (6l[8l. We make the problem 
numerically tractable by combining well-known methods, the Floquet theory for periodic time-dependent 
systems E51 . the Lanczos diagonalization algorithm E6II27 , 28, 29 301, and standard propagation methods 
based on the Runge-Kutta scheme. For both, the diagonalization and the propagation we heavily exploit 
efficient matrix-vector multiplications, storing just the non-zero matrix elements. 

Our results represent benchmarks for the exact numerical treatment of our Wannier-Stark problem and 
similar many-body lattice systems with external drives. This allows for the extrapolation of the applicability 
of our exact methods as a function of memory and computation time, and in view of the physical reality of 
finite size many-body quantum systems. Moreover, our results could be interesting as well for interaction 
effects in the quantum simulation of many particle interference effects with bosonic statistics OT1 . a problem 
known as boson sampling E21 . As a paradigm application of the numerical methods for our Wannier-Stark 
system, we analyse its quasi-energy (or Floquet) spectrum, which is characterized by a large number of 
avoided level crossings in the presence of strong interactions. To do so, we present an efficient and reliable 
way of detecting a large number of such avoided crossings and their distributions, based on the properties 
of the eigenvectors. This allows us to characterize the non-integrability of our model, and, in principle, of 
other many-body lattice systems in general Il33ll34l . 

1.3. Structure of the paper 

In Section[2]we introduce our many-body boson system and we briefly summarize its previously studied 
properties (see refs. H7[|35][36|). The main numerical tools for obtaining the eigenspectrum of our system 
are presented in Section [3] along with a detailed analysis of their efficiency. Section [4] is devoted to the 
quantum spectrum around resonance conditions, for which the two bands are strongly coupled. In particular, 
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Figure 1: (Color online) Sketch of the experimental realization of the system described by (JTJ in a double periodic potential. The 
considered processes of hopping and two-particle interactions are pictured by J a , J h and W x , W a , W/,, respectively, where a indexes 
the lower and b the upper band. The Stark force F, introducing interband transitions, e.g. on-site with a coupling strength Co, acts 
to the left in our sketch. A is the band gap between the two energy ladders corresponding to the two considered energy bands. 


we present an alternative method for detecting and characterizing avoided level crossings. In Section [5] we 
extend the discussion of the spectral properties to the time domain. This turns out to be useful for the 
understanding of the impact of the avoided crossings and their distributions on the system’s dynamics. 
For this, we compare the numerical efficiency of computing the time evolution using either the Floquet 
eigenstates and eigenspectra or a direct Runge-Kutta integration. Finally, Section [6] concludes the paper 
with a brief discussion of the applicability of our methods to similar quantum lattice problems. 

2. Bose-Hubbard model for the Wannier-Stark system 

The Wannier-Stark system is an interesting quantum problem because of its non-intuitive phenomenol¬ 
ogy, including Bloch oscillations (responsible for the effective Stark localization of wave packets). To 
describe interband coupling effects, such as resonant tunneling between different Bloch bands, at least two 
bands must be included into the model. Allowing also for many-body correlations, the problem immedi¬ 
ately becomes very challenging ffT6l [17 ., [9 20,74]. We restrict here to a closed two-band model, which 

(a) could be readily implemented in experiments in good approximation, see the appendix in ifTTl . and 

(b) includes strong interparticle interactions and interband couplings. Our many-body model for bosons is 
based on the celebrated Bose-Hubbard Hamiltonians, usually studied just within a one-band approximation 

uni m m m. 

Our Hamiltonian 


H = 2 -y (p]J l + h. c )+^] 2 P] + ^ Bl hf ) 


IP 


^ cosC^bi + Vyb^aj + h.c. 


If 


£ 21MX+ £(«?-»?), 


( 1 ) 


represents a two-band tight-binding or discrete-lattice model. Its possible experimental realization in a 
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double periodic optical lattice and the considered microscopic processes are sketched in Fig. [T] Our model 
system has the following characteristics: 


(?) it preserves the total number of particles (bosons) N living on L lattice sites or potential wells. Then 
the full dimension of the state (or Fock) space is given by 


, (N + 2L- 1)! 

dp — - 

N\(2L - 1)! 


( 2 ) 


A A+ ( Q\ A, 4- /V 

(?'?) The operators />'/ and /?] represent the annihilation and creation operators, respectively, hj = [3'Jk is 
the number operator, with band index [i e [a, b } for the lower and the upper band, respectively. 

(?'?'?) The parameter space, defined by the parameters {Jp, Wp, W x , C M , A}, is quite large. Jp is the kinetic 
energy scale, characterizing the nearest neighbor hopping matrix elements. The IT’s represent the 
on-site, intra- and inter-band interaction strengths. 2/r + 1 dipole coupling coefficients C /( induced by 
the constant field wj = 2nF, with q e Z. Flence, the Stark force (the tilt) is given by the third term 
in the first line of Eq. Q. The average band distance, the energy band gap, is given by A. All these 
parameters are expressed in the following in units of the recoil energy £ rec = (hk /) 2 /2 m, which a 
photon with momentum hki exchanges when scattering from an atom of mass m. Since all potentials 
are created optically in the experiment, this is the natural scale for realizations with ultracold atoms 
HE!. Furthermore, we set h = 1 throughout this paper. 

(?'v) The available state space T is spanned by the following Fock states 

T = span{|??j ??2 ''') ® \n\n b 2 , •}, (3) 

with YjI n< ? + = N. This space can be split into sets of states with the same upper-band occupation 

number M = (X/ with \d>) an arbitrary Fock state. We refer to these sets of states as M-Manifolds 
in the following. In this way, after relabeling the states using the quantum number M, the Hamiltonian 
matrix takes the following block form 


' H 0 ,o 

H,, 0 

H 2 ,o 

0 

0 


Ho,, 

Hit 





H 0 ,2 




0 

(4) 

0 




Ha/,W-2 





Hat—i,a?_ i 

Ha?,a/-i 


0 


0 

Ha/- 2 , w Ha/-i,a? 

Ha/, n . 



The off-diagonal blocks represent the coupling of the different M-subsets via one- and two-particle 
exchange processes induced by the Hamiltonian terms in the second line of Eq. 0. 

In the non-interacting limit, i.e. Wp_ x = 0, the matrix H can be effectively reduced to a block- 
tridiagonal matrix of size (N + l)x (N + 1) with approximate eigenvalues Em ~ /HA,., see ref. Ii35ll36ll . 
Here 


A r = A ^(1 - cj B r/ A) 2 + 4(m B C 0 /A) 2 
4 


( 5 ) 
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Figure 2: (Color online) Sketch of the structure of eigenvalues for a fixed Stark field F far from the resonance regime, for the 
tilted Hamiltonian in Eq. 0 (lower/blue bars), and of the Floquet quasienergies, i.e. the folded spectrum within a Floquet zone of 
width (x) B (upper/red bars), see sections [2] and |3.1| for details. The values e, are the Lanczos starting points used for improving the 
performance of the diagonalization as presented in section [3~i~| 6s is the maximal splitting of single-particle states generated by the 
interparticle interaction. 


is the energy difference between neighboring manifolds exchanging one particle. The integer r labels 
the order of the resonances lUTIl occurring for tilts F r ~ A/2 nr. These resonances appear at certain 
values of the force for which the single-particle levels in the upper band become degenerate with 
the ones of the lower band dim The region in parameter space as a function of F close to F r 
is denoted as resonant tunneling (RET) regime in the following. The interactions, together with 
the strong interband coupling at RET, lead to a strong mixing of levels. In contrast, we can resort 
to a perturbative argument far off the resonances. Off-resonance, the eigenvalues Em are corrected 
by including the splitting of the single particle levels due to the interactions. The corresponding 
eigenstates can be labeled by the integers 0p=\ t 2 — (3 1 ))<p and 6 X = 2(£/ Therefore, 

they can be approached by the formula IIT71 [351 : 


% = MA r + G ■ W, 


( 6 ) 


with W 7 = {W a , Wb, W x ), and 6 = (0 a ,0 b , 0 X ). The maximal intramanifold splitting generated by the 
interaction can be estimated by 8s = rnaxj U^ h }, with 


(U a )max 

w 

= -^-(N - M)(N - M - 1), 

max 

W b 

= yM(tf-l), 

(U%) max 

= 2 W X (N - M)M. 


The values of Uff b correspond to the largest energy cost for producing states with M particles in a 
single lattice site, e.g., 10) ~ \N - M 00, •••)<£> \M00, ■ ■ •). In consequence, in most off-resonant cases 
with F far from F r , we have A,- ~ A » W a j KX , and, therefore, the overlapping between contiguous 
M-manifolds can be neglected. Then, the spectrum consists of bunches of levels within an energy 
range 8s around Em, as sketched in Fig. [2] 

According to Eq. ([ 6 ]), the intramanifold energies Em,b, far from the F r , are proportional to the tilt, c.f. 
also Fig. [3j a). Therefore, the parametric dependence of the levels on F, s(F), is essentially given by 
straight lines parallel to each other, with slopes w -MA/F r . Conversely, the intermanifold energy 
levels approach each other in the vicinity of the interband resonances at F = F r (see Fig. [3]). Here 
avoided crossings take place involving all the levels arising from the lack of additional symmetries 
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and from the strong interband couplings ifTTl . Within the resonant regime, Eq. ([ 6 ]) is clearly no longer 
valid. The eigenenergies lose their manifold structure since typically we have A » W a ,b,x ~ A 
Nevertheless, Eq. ([ 6 ]) is still useful for performance tests of the numerical routines presented in the 
next sections. 

(v) A helpful trick for the numerical solution of the Hamiltonian (fT| is to transform it into the interaction 
picture with respect to the tilting term • This procedure allows us to recover the transla¬ 

tional invariance in the new Hamiltonian H'. We can, therefore, impose periodic boundary conditions 
in space, i.e., fi L+l = j. Thereby, we now work with the common eigenbasis of the translation op¬ 
erator 5 and the transformed Hamiltonian H' since [ H ', S | = 0. This new basis is the translationally 
invariant Fock (TIF) basis indexed by a and introduced, e.g., in ref. lfl9ll : 

D a 

\s a ,«) = D~ m 2 e i2ml S l \ny‘ • ■ ■> <g> | n\n\, ■■■>«, (8) 

/= t 

where S m |s a , k) - e~ l0JBKm \s a , k). The action of the translational operator on a Fock state is 

s m \ni •••«/•••> = | m +m ■ ■ ■ m+ m ■ ■ ■). 


k = Kj = j / D a ( j 6 |0, D ir - 1]) is the lattice quasimomentum and D„ is the total number of cyclical 
permutations of the Fock state. The dimension of the Hilbert space expanded by the TIF basis is 
d s ~ df lL. It is easily seen that the block structure of the Hamiltonian matrix remains invariant and 
the M- sets, now j |,v f> ., k, M)}, have dimensions 


1 (M + L - 1 
dM = L\ L - 1 


r v.r\ 


(9) 


where (• • •) stands for the combinatorial function. The full dimension of the Hilbert space is then 
given by d s = Zm d m■ 

A second important consequence of transforming H into H' is that the new Hamiltonian is explicitly 
time dependent. This originates in the transformed terms 

-» P]+\Pi exp (-ioJrit) 

ci] + b, -> d] + J>i exp (-icoBUt). ( 10 ) 

One can easily see that H'(t) is periodic in time with a fundamental period Tg = 2k/ojb. the Bloch 
period. The temporal periodicity allows one to obtain a stationary spectrum, not of the Hamiltonian, 
but of the single-cycle Floquet operator llT9ll25 1 


U(T b ) = f'exp|-i J B dtH'(t) 


( 11 ) 


where T implies time ordering. The eigenphases cp j 


U(T B )\<f>j) = , 


( 12 ) 


divided by the period T B define the quasienergies, in analogy to quasimomentum for the Bloch prob¬ 
lem of a particle in a periodic potential. The set of quasienergies {e,} inherits the periodicity a> B in 
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Figure 3: (Color online) Wannier-Stark spectrum for the single-particle case (a,b) and typical system parameters. Panels (c) and 
(d) show the resonant regime for N/L = 4/5 and weak (W a = W h = 1.3 x 10 -3 , W x = 1.4 x 10~ 3 ) or strong ( W a = Wj, = 
1.3 x 10 2 , W x = 1.4 x 1(T 2 ) interparticle interactions, around the second order resonance at tn s /A = 2 for J a = 0.059,/* = 
-0.071, Co = -0.095, Cj = 0.043, C 2 = -0.005, A = 0.56. The color code is as follows: eigenstates with M t ^ 0 (black lines), 
eigenstates with M, = N (red/dark grey lines), and eigenstates with 0 < Mj < N (green/light grey lines). For better visibility of 
the avoided crossings, we plot the panels (b-d) directly as a function of the Stark force F oc u) B , whilst the fan of single-particle 
Wannier-Stark levels in (a) is best seen as a function of the inverse Stark force oc l/w B . 


energy space. This is a direct consequence of the Floquet theorem for periodically time-dependent 
systems fl25l (291 EL Therefore, we can choose, without loss of generality, ej e (—a> B /2, a> B /2], 
implying that the spectrum is folded within an energy range of size wj, the fundamental Floquet zone 
(FZ), as sketched by the red bars in Fig. [2] 

The so-called extended spectrum is obtained by shifting the eigenenergies from FZ as Sj —> sj+rifiOB, 
with the integer iif labeling the FZ. To obtain numerically the spectrum directly from U(T B ) typically 
demands long computation times since the often applied approach consists of two steps. First, one 
computes the matrix U{T B ) by propagating the entire basis, for instance, the TIF basis, from t = 0 
to t = T b . Subsequently, the obtained matrix must still be diagonalized, whereby it is not a sparse 
but a full matrix in general | |T9l TOl . In this paper, we follow a different but equivalent approach: We 
study and diagonalize the Floquet Hamiltonian to be introduced in the next section. This procedure 
avoids the time integration for obtaining the one-cycle operator from Eq. as done in previous 
single-band studies of many-body Wannier-Stark problems Ifl5l[l9ll20] l. Now we directly diagonalize 
much larger, but sparse matrices instead. We shall use the properties of our Hamiltonian 0 described 
above in order to improve the performance of the diagonalization scheme that will be presented in the 
next section. 


3. Floquet-Lanczos Diagonalization: Generalities and performance 

In this section we discuss the implementation of two well-known techniques, Floquet theory and Lanc- 
zos diagonalization. The combination of these two methods allows for a detailed numerical treatment of the 
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Hamiltonian H' describing our many-body Wannier-Stark problem for reasonable system sizes. 

3.1. Floquet theory 

Let us start with the Floquet representation of transformed H'{t). The periodicity in time of H\t), i.e., 
H'(t + Tb) = //'(f), permits the application of the Floquet theorem |[25l [29l [37l . The Floquet Hamiltonian, 
defined as Hf = //'(f) - id,, has instantaneous eigenvectors 


Hf\sj{t)) = e j \s j {t)), (13) 

which have the same periodicity as //'(f), i.e. |e,(f + Tb)) = |e,-(f)). Any time-dependent solution can be 
expressed now as 

tyt) = ^ exp (-isjtj | Sj(t )), (14) 

j 

where the coefficients Cj are the projections of the initial state onto the Floquet basis vectors. 

In order to obtain a system of stationary equations, which we can solve by diagonalization, we expand 
these vectors into their Fourier series (making perfect sense because of their periodicity). The corresponding 
multi-mode Fourier decomposition reads 

\ej(t)) = J] e~ Aat \4j) = 2 e- lk01,< %) • (15) 

i k 

The expression on the very right is possible because of the existence of a fundamental frequency cob since 
also the remaining frequencies are multiples of w#; so we can write: k ■ a> = (k\ + 2k^ + ... + = ka)/;. 

The eigensystem of equations is then given by the following coupled algebraic equations: 

eA\4j) = (A) - ojBkl) 14 ) + j\<p k ~ l ) + 

+ ^ [z,^) + z;\4^)}. (i6) 

Here, Hq collects all the time-independent parts of H\t). The operators J and Z /( define the hopping and 
dipole-like transition operators 


1 V 1 ^4- 

1 = -nYjWi 

IP 

Zu = ^sC^^a^bi. (17) 

i 

Equation ([16]) can be rewritten in the matrix form 

(M - SjS) ■ xj = 0, (18) 

after expanding the states of the kth Fourier component | <p*.) in the TIF basis, |0^> = Zio- C J ka \s a ), and 

defining the vector of components x r . = (c J k a , ■ ■ ■ , C J k a , ■ ■ ■, C J k j. In our case, 23 = 1, giving a standard 
eigenvalue problem. The resulting square matrix is banded and blocked and it has the total dimension 


«tot = d s Ak , 

8 


( 19 ) 


where A k = k max + k min + 1 is the total number of Fourier components considered, with -k, mn < k < k max . 
The matrix At is thus given by 


M = 


H 0 + ^max (Ob 1 


H 0 + 2mgl 
J + Z, 

Z 2 

Z 3 


J f + z{ 

Ho + o»b1 

J + Z! 

z 2 


A 

J f + z\ 

H 0 

J + Z! 


A 

A 


J f + z\ 

Ho - tufil 

J+Z! 


3 

z t 


J + z\ 

Ho - 2o>b 1 


Ho 

( 20 ) 

The matrix Al inherits the sparseness of H; the number of non-diagonal blocks depends on the number of 
considered //-dipole couplings. In the current case, we can restrict to |//| < 2, for which // = {—2, —1,0,1,2} 
since the dipole coefficients go fast to zero with increasing |//| lf36l . As a consequence, the matrix At 
is effectively a very sparse, five-diagonal block matrix. As for the operator U(Tb), the eigenvalues of 
the Floquet matrix lie within the fundamental Floquet zone sj e (-cob/ 2,o»a/2], with £/mod(a/g). The 
eigenvectors are obtained as the linear combination | Sj) = TjkaC J ka \s a ) = C J a \s a ), with = 0*- 

The dimension of the energy basis of interest, {\£j)}j, is therefore d s . 

Since a maximal number d s (<sc /i to t) of eigenenergies lie within one FZ, the computation of the full set of 
eigenvalues of At («tot) is not necessary. Therefore, it is convenient to implement a sophisticated algorithm 
that permits us to compute a controllable number of eigenvalues and eigenstates of At. We implemented 
the symmetric Lanczos algorithm, which proved to be very efficient for obtaining either all or at least a 
reasonable number of eigenvalues in the fundamental FZ. Without loss of generality, we work with the TIF 
subspace defined by quasimomentum k = 0 (see Sec. [2]). This choice automatically guarantees Al to be a 
real and symmetric matrix. 


3.2. Lanczos diagonalization 

The sparseness of the matrix At makes it very suitable for numerical diagonalization by means of 
algorithms such as the symmetric Lanczos procedure ESI . This method allows us the computation of 
eigenvalues and eigenstates of a given complex symmetric matrix. A nice feature of this method is that 
it permits the computation of eigenvalues close to a given initial choice in energy space. Let us start by 
shifting the diagonal elements as Al —» Al + AfL The goal is to find the largest eigenvalues Aj = 1 /(sj - £o) 
of the matrix At -1 , which correspond to the eigenvalues sj close to the predefined shift parameter eq. 

The method is an adaptation of the power method |[28ll to find eigenvalues and eigenvectors of a square 
matrix based on the construction of the Krylov subspace (A„ ;) ); that is, the space spanned by 


%i p {M 1 ,77i] = span{7 /!,/7 2 ,---,?7 „ p }, (21) 

with T]j+ j = 7/i is a suitable initial vector, which is normalized ||t/i|| = 1 . n p («: n t ot) is the 

dimension of the Krylov subspace, here referred to as Lanczos iterations number. We define the matrix 
P ~ [ 7 / 1 , 77 ?,..., 7/„ ;j _i j. The Krylov subspace representation of At -1 is the matrix T = 'P ~ { which is 
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tridiagonal and symmetric, with elements given by 


a i ■■ T ij ~ 'IjM '///. 

fij T jJ+ 1 = ||tj +1 t /+ i||, (22) 

with 

tj +1 = M ~ 1 rjj - ajT]j - (Sj- 1 Tjj- 1 

iJj+\ - tj+i/Pj+i- (23) 


The construction of the < 7C„ space is supported by an appropriate Gram-Schmidt orthogonalization proce¬ 
dure for the rn vectors. The matrix T can be diagonalized using standard diagonalization routines. Here we 
use the QR algorithm Il28l with an error of the order of the numerical (double) precision I (T 14 ; that is, we 
compute the decomposition T - QR and the subsequent iterations 

r M = RjQj = Q’jTjQj = QfTjQj , (24) 


with To = T. The matrix T j+ 1 , by the equation above similar to T, converges to the Schur form of 
symmetric matrix T, i.e., the eigendecomposition of T. In addition, the discrimination criterion introduced 
by Parlett and Scott E71 is implemented in order to separate converged from non-converged eigenvalues. 
Finally, since the matrices T and At 1 are similar, they have the same eigensystem. 

The numerical implementation of the Lanczos algorithm requires the setting of several parameters. The 
first of these is the shift so, which we set using Eq. for instance, so - Em- This allows us to find a 
set of eigenvalues in the vicinity of manifold energy Em = MA r of interest to us (see Fig. [2]). Using the 
approximately known dependence of so(F ) on the Stark force F permits us to set the center of the Floquet 
zone, for instance, in middle of the spectral range, so ~ j/VA, for convenience. This choice helps in the 
visualization of the spectrum in the vicinity of F r , as shown for instance in Fig. [3] 

The number of the found eigenvalues depends on the number of Fourier components A k and the number 
of Fanczos iterations n p . k max can be estimated by the ratio of the eigenenergy bandwidth of H in Eq. 0 
and the size of the FZ o>b- In this way, k max can be regarded as the number of Stark-induced ’’fictitious” 
photonic excitations Hcob necessary to promote all atoms from the lowest (M = 0) to the largest (M = N ) 
manifold excitation (see the sketch in Fig. [2]). This value is thus given by 


k 


max 


A E 

U>B 


NA r + W h N 2 
2:tF 


(25) 


where is the energetic cost of producing a state with N atoms in a single lattice site within the second 
band |[T7ll . e.g., the state \sj) ~ |000- • •) 0 |/V 00 ■ • •). For our purpose, it is very important to note that 
according to the estimation of k max , we have that k mui does just play a minor role. Therefore, the Fourier 
expansion may not necessarily be symmetric with respect to the photon index k, i.e. k n \ m /k max < 1 is 
indeed a possible choice for us. To minimize the dimension of At, we may even choose k min = 1 and 
still obtain converged Floquet eigenspectra. The number of Fanczos iterations n p is adapted such that we 
obtain a minimum number d s of convergent eigenvalues within the FZ, for which we usually have that 
d s < ll p « 72tot ■ 

The required memory storage for one diagonalization is given by the formula 

MS[GB] = (/2tot72iarg + «£) X ~ d] , (26) 
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Figure 4: (Color online) Spectra computed at the first order resonance between the two energy bands for increasing sizes of 
our system, i.e. for the filling factors shown in the legend in (b). Panel (a) shows the saturation of the number of convergent 
eigenenergies (#CEvs), within one Floquet zone (FZ), with increasing size of the Krylov subspace. The y-axis is normalized to 
the corresponding TIF basis dimension, d s and fc max , are estimated as in Eq. m- Panel (b) shows the saturation of #CEvs with 
the number of Fourier components A k. The realistic parameters IH][36) for this computation are: r = 1, A = 2.53, J a = 0.08, 
J b = -0.24, W a = 0.018, W b = 0.025, W x = 0.02, C 0 = -0.09, C ±l = 0.035, and C ±2 = 0.002. 


where / 2 | arg ~ (ju + 1 )d s and n p = dim(7C, ,)■ The first term on the right hand side is the needed memory 
storage for At, which for systems with fillings, i.e. particles per number of sites, N/L > 1 can be very 
large. n p is the dimension of the Krylov matrix P, which is proportional to the dimension of the TIF basis. 
For storing a 64bit number we need 8 bytes. The factor 2x8 = 16 comes from the fact that we use an 
implementation of the Lanczos algorithm which was developed for general complex symmetric matrices 
| |29j , [30] [37l l . For instance, for N/L = 7/6 (d s - 5304) we need k mdX > N/r Fourier components at RET 
conditions F = F r , given A » ( Wp tX , C /( |. Therefore, we have in this case that MS[GB] ~ 35. The computer 
used for these calculations was a 24 core AMD Opteron(TM) Processor 6234 with 60-GB main memory. 


In Fig. §(a), we show the dependence of the number of convergent eigenenergies of A4, normalized 
to d s , on the size of the Krylov subspace, for different filling factors N/L. Once k max is set by Eq. (25), 
we see that the saturation point ( n p /d s ~ 1.3) is the same for all systems. In Fig. j4j-(b), we scan instead 
the number of Fourier components by varying k max for fixed size of the t K„ p subspace n p = 2d s . Figure [5] 
shows a polynomial growth of the computational time (CPU time) with increasing dimension of the TIF 
basis. We observe a polynomial increase of the CPU consistent with ~ d 3 . Indeed the only limitation in our 
numerical scheme is the memory storage since k max ~ u> B l = \/2nF, which implies that the dimension of 
the At increases considerably for smaller and smaller values of F. This fact makes it difficult to study, for 
instance, high order resonances with large r, for which the tilt is small. A positive aspect of our system is, 
however, that many of its properties are not really size-dependent, but they rather depend only on the ratio 
N/L (see ref. ifTTI ). The growth law for the computational time and, in consequence, the efficiency of our 
numerical method are expected to hold for even large systems as may be extrapolated from Fig. [5] 

One way to overcome the memory storage limit, e.g. for small values of F, is to simply compute only 
parts of the spectrum in single runs. As we can see from Eqs. ([5]) and ( |25| ), and also from the manifold 
structure of the eigenspectrum folded within the FZ (see sketch in Fig. [2]), we can improve the performance 
by setting various starting parameters Sj for the Fanczos algorithm, while defining k nVdX = k nun ~ A r /2o»g. 
We now can proceed by implementing N + 1 independent Fanczos diagonalizations, each with its own 
sufficient number of Fourier components to promote one atom to an energy difference A r /2 starting from 
every Sj. This only involves excitations until half way between two neighboring manifolds M and M + 1. 
Thereby, we expect to obtain approximately a number ~ dM of eigenenergies around every sj. In this 
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Figure 5: (Color online) Log-log plot of the computational time (CPU time) vs. the TIF basis dimension for the quasimomentum 
block with k = 0. We extract the following relation for the CPU time ~ d" with n = 3.2 ± 0.2. The exponent n is estimated by the 
slope of the fitting function for the data corresponding to N = 5 (black o), given by the black straight line. Note that for N = (6, 7) 
no relevant changes occur as one can see even by eye. d s is increased by fixing the total particle number N and varying the number 
of lattice sites over one order of magnitude L = (3, • • ■, 12). The dimensions of the corresponding Floquet matrices At are then 
n, ot — 10 2 - - ■ 10 5 . The system parameters are the same as applied in Fig. [ 4 ] 


scheme, the dimension of the Krylov subspace ( n p < d s ) is much reduced along with the dimension of 
the M. Consequently also the memory storage and the computational times are reduced for the individual 
runs. In order to obtain the full spectrum around any chosen shift, e.g., £0 ~ 0 as in Fig. [2] we wrap all 
eigenvalues (obtained around the others shifts ~ 0) into the FZ defined by so as 


Si - koa>B , with ko - 


s j± 0 + d — so 


Cl>b 


eZ- 


(27) 


Here [■ • •] stands for the integer part function and ko is the number of Floquet zones needed to map e'.s onto 
their corresponding Floquet eigenstates within the fundamental FZ centered at so ~ 0. We defined 6 = 5s 12. 
The folded eigenvalues thus satisfy sj = epnodfco/}). As the final step we discriminate and eliminate copies 
of eigenstates using the fact that the overlap of physically different eigenvectors ideally vanishes. Our 
numerical orthogonalization criterion is that if Ke;|e/)| < 10 “ 7 we assume that the two vectors |e,) and | sj) 
(i + j) represent different eigenstates. 


Table 1: Comparison of the eigenenergies of the full spectrum (*) around So ~ 0 (upper table) and e 2 ~ 2A r (lower table) with 


folded eigenenergies with starting Lanczos parameters Sj = j A,- (+c>) (*'°), with 8 as defined after Eq. l27 1 . The Stark force is 


F = Ax 10 3 In, while the other parameters are the same as used in Fig. [ 4 ] 


£o 

0.001349897362263* 

0.000531659341077* 

-0.023753576110081* 

-0.011655535562532* 


£ 1 *, £1 + 6° 

£ 2 *, £2 + d* 

£ 3 *, £3 + d* 

£ 4 *, £4 + d* 


0.001349897362262* 

0.000531659341080* 

-0.023753576110153* 

-0.011655535570281* 


0.001349897362263* 

0.000531659341075* 

-0.023753576110091* 

-0.011655535561856* 

£2 

5.040989227022449* 

5.083166308393272* 

5.071405551272533* 

5.056877702765119* 


£ 0 *, so + 

£ 1 *, £ 1 + d* 

£ 3 *, £3 + d* 

£ 4 *, £4 + d* 


5.040989227022449* 

5.083166308393297* 

5.071405551272531* 

5.056877702765120* 


5.040989227022451* 

5.083166308392099* 

5.071405551272880* 

5.056877702764962* 


12 


















In Table [T] we present an example of results based on the implementation with varying Lanczos shift 
parameters. We computed the full spectrum around so ~ 0 for N/L = 5/5 with d s = 402, see the given 
eigenvalues marked by (*). For this calculation, we used the parameters A, * A = 2.53, 2nF = 2A x 10 -3 , 
and Wb = 0.025. Therefore, the full diagonalization requires k max « (N/2 + 0.125) x 101 « 262 (with 
fixed k m i n = 1) Fourier components to obtain d s converged eigenvalues. The computation time for this 
computation was CPU time w 20 minutes. In comparison, we computed the subsets of eigenvalues around 
the shift parameters e ; >o = {A r , • ■ • ,NA r }, with k max = k rnin ~ A/2 lor = 101/4 w 25 Fourier components. 
Then, by using Eq. (27 1 , we folded them into the FZ centered at eq ~ 0. For a single diagonalization process, 
the CPU time ~ 95 seconds, that is, the total CPU time is « 95 x (N + 1) seconds « 10 minutes. This is 
at least by a factor of two faster than the time needed for one full diagonalization of the same problem. 
Hence, by using different shift parameters, a clear reduction of the computational time is expected for large 
system sizes. Then even computations for small tilts F become possible. The resulting folded eigenvalues 
are shown in the remaining rows of Table [I] (eigenvalues (*’°)) in comparison with the ones from the full 
spectrum around eq « 0 (*). Note that the agreement is excellent up to 12 digits, at least for the first three 
eigenvalues belonging to the shift parameters {e 7 =i, 2,3 = ./A, (+d)}. There is a growing lack of accuracy as 
the energy distance increases between the FZ of interest (eo) and the one used for calculation of the spectrum 
increases (sj± o). This problem can be partially corrected by choosing a different starting Lanczos parameter 
for the FZ of interest. That is, we redefine eq £2 = 2A r , and again compute the eigenvalues around 
{is# 2 }, and fold them using Eq. (27 1 . The result is shown in the lower part of the table. Here the agreement 
even with the largest eigenvalues, corresponding to those of the high lying manifolds, is much improved. 
For our specific Wannier-Stark system, the trick of using different shift parameters is well supported by the 
quasi symmetry related to the occupation number M. This symmetry is based on the fact that, for F = 0 
and small W x compared to the manifold band gap A,, the Hamiltonian can be approximated by the block 
form H = Hm. As explained in Sec. [2J this quasi symmetry remains approximately valid for F + 0, 
but only far from the resonance conditions F - F r . 

So far we have presented an efficient scheme for the numerical diagonalization of the Floquet Hamilto¬ 
nian represented by the matrix At. With both eigenvalues and eigenvectors at hand we shall now study the 
structure of the spectra with respect to changes of a control parameter, which is naturally the Stark force F 
in our case. The extension and connection of these static properties of the spectra to the temporal evolution 
will then be discussed in section[5] 


4. Detecting avoided crossings by eigenstate projection 

After solving the Floquet eigensystem with the methods exposed in section [3] we now analyse the spec¬ 
trum of our non-integrable many-body system. The complexity in the spectra arises from the level-level 
repulsion of eigenstates as a function of, e.g., the control parameter F. Strong couplings in the resonant 
regime, in the vicinity of F r ~ A/2/zr, typically lead to avoided crossings (ACs) in the parametric depen¬ 
dence of the levels. The occurrence of ACs with a wide distribution of crossing widths is associated with 
the onset of chaos in quantum systems (see e.g. the refs. ifTTl [T9l l35l l38ll in this context). The chaoticity of 
our system can be probed by means of statistical distributions, for instance, the nearest neighboring spac¬ 
ing distribution Il38ll39l at fixed values of F. Another way to do this is by studying the distribution of ACs 
widths within a small range A F Ii40ll41 1. In the following, we focus on the latter approach as it characterizes 
a paradigm phenomenon which is directly connected to the eigenstates of our system. 

The locations of ACs and their widths are difficult to detect numerically in large systems. Most methods 
are based on following the parametric changes of the curvature of the energy levels If34l . or on solving 
non-linear equations for the discriminant of the eigenvalue equation for a given Hamiltonian matrix lf42ll . 
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The implementation of these methods affords the ordering of the energies and, for instance, the implemen¬ 
tation of the discriminant method which is computationally expensive for large system sizes. Here, we 
introduce an alternative method for detecting and characterizing ACs based on the statistical properties of 
the eigenvectors. The method is suitable for our Floquet system with periodic quasienergies since it avoids 
the implementation of routines for sorting and following the quasienergies. Our procedure may be applied 
to any kind of discrete quantum spectra making it an interesting method for the spectral analysis of generic 
quantum systems. 

Relevant information on a complex quantum system can be extracted from the properties of the eigen¬ 
vectors of the Hamiltonian under consideration Il43l 1441 . Let us illustrate the idea with a simple example 
in our case. The eigenstates | sj) are connected to the TIF basis through a unitary transformation matrix 
as e = 'ZFs, with e 7 =(•••, | sj), ■ • •) and s T = (• • •, • • •)• The transformation matrix elements are 

(^)aj = (sa\£j)- In this way, since s = F/e, any TIF vector can be written as |.v„> = Yij( s j\ s a)\£i)- 
Since the occurrence of ACs is typically local in the spectrum, there are only a few eigenstates involved, 
two in the simplest case. Therefore, we may approximate ~ (e ; |5 Q: >|£j) + (gy+i ), where 

Ke/K)! 2 ~ Key+ilsQ,)! 2 . If now Ke / i.s' a >| 2 + Key+iUo.)! 2 ~ 1 is satisfied, the state |s a ) approaches the hy¬ 
bridized state 

l4 ±) > ~ 4= (ley) ± |e/+t>) ■ (28) 


Such a hybridization in the diabatic basis is typical of an AC. It effectively means that the effective number 
of participating states in the representation of the local eigenstate is larger than one. Quantitatively, we can 
express this using the inverse participation ratio |j43l 

ipr a (F) = ^ l<£,K>| 4 , (29) 


which is then also larger than one. Note that, for the state in Eq. (28 1 , we have l/ipr a = 2. Sufficiently far 
from an AC, ipr ~ 1 < 2, which implies that 1 /ipr ff has a local maximum at the exact position of the AC. The 
width of the AC is the energy separation c = s- l+ \ - £,. Note that this simple example does not involve the 

Ql 2 . 


eigenenergies at all, but it does need the projections p a . = \(s 
c J 


j\*a/ 


The method can be generalized as follows. Let [)J y be the vector of projections p J(F) = (■•■, p°j{F), ■ ■ -j. 
It represents the TIF state |.v„) in the energy basis. Given two of these vectors, p and q, they are never 
collinear unless they are the same, i.e., p = q. Therefore, any change in the structure of p„, as a function 
of the control parameter F, corresponds to the coupling to additional states. This is exactly what happens 
at ACs. To quantify the structural transformation of p, T , a well suited measure is the fidelity function 
introduced in OT 1 as 


G a (F, 6F) 


pI(F)PJF + SF ) 

pJ(F)|l+dE^Jp a (F). 


(30) 


Here we have used a Taylor expansion up to the first order in 6F\ 

P a(F + SF) = p a (F) + SF—p a (F ) + • • • . (31) 

oF 

The differential operator Tgp = 1 + SFd/dF acts as a translation operator in the F-coordinate. It obeys 
the composition property, i.e., Tsp l Tgf 2 = T l} r l+ sri, with Tq = 1. Using the orthonormality condition, 
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Figure 6: (Color online) Single-particle spectrum (a) and the strong coupling at the ACs visible in the effective number of partic¬ 
ipating states 1/ipr (b). The color discrimination is as follows: eigenstates with M,- - 0 (black lines) and eigenstates with M, =; 1 
(red/grey lines). The band gap is A = 2.53. 


p'j = 1, it is straightforwardly shown that 

llpj := a/p^P o' = y/d s cr 2 + 1 /d s , (32) 

with cr 2 being the variance of p ff . Rewriting the fidelity function from Eq.(|37]), we obtain 

d, 

G a {F,6F) = T sm Yp'p 2 ' (33) 

7=1 

Here the last term, on the right-hand side is nothing but the inverse participation ratio ipr ff (F) = 
Z j |<$ ff |e/F)>| 4 . Its inverse defines the effective number of participating states, just as in the example above 
for two states. Finally, we arrive at the relation 

lim G a (F, 8F ) = ipr ff (F) = d s cr + d;' . (34) 

<SF—>0 

The latter implies that we have to find the maxima of the 1/ipr in order to locate the position of the 
AC. The upper and lower bounds are 1/ipr = 1 and l/ipr ff = d s corresponding to the vectors p ( ^ = 
(0,0, • • •, 1, • • • ,0,0) and = d; 1 (1,1, 1, 1,1), respectively. This can be seen in Fig. 6, where 

the lower panel shows values between 1 and 2 for ACs with two coupled states, c.f. also the description 
after Eq. (29). 

To estimate the width, c, of the detected AC, we appeal to the fact that the occurrence of an AC is 
typically local in energy space. Thereby, it is expected that at least two neighboring components of p tf are 
maximal whenever the hybridization process happens. We can then study the functional behavior of the 
local energy difference co = cj = sj+ \ - sj by means of the following projection vector distribution for two 
levels 


Pa\od) = J] P^P^+l • ~ ( £ 7+l “ S J 0) 

j 


r-o n (a> - (e 7+ i - sj )) 2 + T 2 


( 35 ) 
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Figure 7: (Color online) (a) Many-body spectrum around the first-order resonance, r = 1, for N/L = 4/5 and the k = 0 TIF 
subspace. The zoomed region on the right highlights the clustering of avoided crossings in the center of the spectrum. The 
density of ACs is shown in (b). Panel (c) presents the normalized distribution of the effective number of participating states, 1/ipr, 
computed at all detected ACs within A F = 0.5 x A/2n, i.e. in the full range of F shown in (a). All panels (a-c) visualize the strong 
coupling, most prominent at the center of the spectrum. 

The system parameters motivated in 1361 are: A = 0.56, J a = 0.06, J/, = -0.071, W a = 0.025, W/, = 0.027, 

W x = 0.024, C 0 = -0.095, C ±1 = 0.04, C ±2 - 0.004. 


(2) 

Here we expressed Dirac’s delta in terms of a Lorentzian function with vanishing width T. P a (a>) as a 
single function efficiently detects the widths a> = sj+i - sj at which P„. (o>) has a maximum. In a more 
general case, an AC can consist of more than just two participating states. In this case, for a better estimation 
of the AC positions, we also use the following three level distribution P^Hoj) = X, //*_ ] p"p'j + ] 6(u - c), with 
c = minjey+i - Ej, ej - Ej. |}. Because of the typically local feature of ACs in a discrete quantum spectrum 
(see the zoom in Fig. 17] as an example), it is practically sufficient to rely on the information content of 
both functions P a and P a for a reliable detection of ACs. We used this approach to find numerically the 
ACs, their positions, and their distributions within the resonant (RET) regime in an efficient and easy to 
automatize manner. 

In figure |7] we present the results for the detection ACs in our many-body Wannier-Stark system. In 
the zoomed region one can recognize the strong level repulsion between the energy levels as a function of 
the force sj(F). The ACs cluster in the region around F r , as highlighted by the density of ACs Pac(F) - 
#ACs/d s AF in Fig. ]7jb). Furthermore, we show the normalized distribution of the effective number of 
participating states in the entire resonance region A F in Fig. ]7Jc). We see the large number of effectively 
participating states in the resonant regime, arising from the strong level mixing, in particular for the states 
with manifold number 0 < M < N (see section |2]). 

As shown in ref. ifTT! for our system, the number of ACs can be increased in our model by increasing the 
filling factor ( N/L >1) and the interparticle interactions. This manifests itself in a compression of levels 
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Figure 8: (Color online) The inset shows the comparison between the theoretical predictions based on random matrices for P(c) 
(red/dark grey dashed line) and the numerically obtained distribution (black line). The main panel shows the respective cumulative 
distribution functions. Good agreement between the data and the random matrix theory for quantum chaos, expressed in the 
spectral statistics, clearly highlights the strong coupling and, consequently, the complex dynamics in the temporal evolution of 
the system. Here N/L = 5/3 (d s = 84) and N/L = 7/4 ( d s = 858) with a bandgap A = 0.155 guaranteeing a large number of 
ACs (#ACs = 4751 for N/L = 5/3 and #ACs = 152145 for N/L = 7/4). AF is chosen as for the previous figure and a step size 
6F = 1CT 4 was used for the computations. 


within the Floquet zone ifTTl . Due to the lack of additional symmetries in the system, this leads to a crossover 
between regular and quantum chaotic spectra and dynamics, where the chaotic regime occurs for N/L ~ 1 
and A < 1 in the vicinity of the resonances at F = F r . To characterize the quantum chaoticity, we need not 
only the pure relative occurrence frequencies of ACs, but also the statistics of their width distribution. We 
recall that our Hamiltonian in the k = 0 TIF subspace is represented by a real and symmetric Floquet matrix, 
which, within the context of random matrix theory (RMT), belongs to the so called circular orthogonal 
ensemble (COE) If38l . Random matrix theory predicts for this case the following distribution of AC widths 
form fiTl 


2 y 2 I y 2 c 2 

P(c ) = (1 - r)S(c - Co) + — exp 

n \ n 


(36) 


with normalized average distance (c) = 1. y is a fitting parameter. For y = 0, the formula characterizes 
a fully regular spectrum, with a relevant energy scale co- Conversely, y = 1 implies the randomization of 
the set of widths c, marking the completely chaotic regime (see Ref. lf4T1 for details). Here we confirm the 
expected chaoticity of our system by comparing our numerical distribution with the theoretical prediction 
for N/L = 5/3 (green/light grey) and N/L = 7/4 (red/dark grey) in Fig. [8] Interestingly, we see that even 
for a small system N/L = 5/3 (d s = 84), the agreement between the chaotic distribution with y = 1 and 
the numerical one is still good. The agreement is excellent for the larger fillings N/L = 7/4 (d s = 858), 
which guarantee better statistics. The number of detected ACs is 4751 for N/L = 5/3 and is 152145 for 
N/L = 7/4. The quality of the matching of the results of our numerical experiments with the theoretical 
predictions is best corroborated by computing the cumulative distributions 1(c) = f Q P(c')dc ', which are 
shown in the main panel of Fig. [ 8 ] 

To summarize, we have shown that the detection and characterization of ACs is possible using the 
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projection vectors p„ and their transformation as a function of a control parameter. The presented numerical 
method is well suited for any class of discrete eigenspectra, that is, not only for our Floquet spectra. Since 
the ACs occurrence is mostly local in the spectrum, the use of the two and three-neighbor level distributions 
turned out to be sufficient for well estimating the width of the ACs, such as shown in Fig. [ 8 ] The parametric 
evolution of the projection vectors studied here is inherently connected with the dynamical evolution of the 
system. This will be seen in the next section after we will have discussed numerical techniques to evolve 
arbitrary initial states in time. 

5. Temporal evolution: algorithms and results 

Flaving at hand the solution of the eigenvalue problem for our Floquet Flamiltonian in Eq. ( fi~ 8 j ), we can 
easily propagate any given initial state \ifjQ) in time. To do this, we best remind the reader on the representa¬ 
tion of the time-evolution operator in the eigenbasis. With help of the closure relation 2/ |£/(0)( £ /(0l = 1, 
the evolution from time t\ to time tj is obtained by applying the following operator onto the initial vector 

d s 

U(t2,h ) = J]e-^ Af |£/f 2 )><£;(?l)[ 
j= 1 

ds ^rnax 

= Z Z (37) 

7=1 kk f =—k m i n 

where At = U - t\ f29l . We can then compute any relevant physical quantity using this operator and the 
coefficients c 9 = Yiki&Sjl'I'o) defining the initial condition. Note that, in the case of = |s a ), we have 

c j = £m( c L)*Wo) = Yj (c L)* = • os) 

ka' k 

This straightforwardly shows us the direct connection between the initial condition for the time evolution 
and the projection vectors defined in the previous section. Starting the evolution at t = 0, we can compute 
any observable described by a hermitian operator O using U(L 0 ) as follows 

o(o = mm) = 2( a °(0)*aJ(o<^iom , m 

ap 


with 


A «(0 = YjCy i{e ’ +0JBk)l CL- (40) 

jk 

Since for large systems, with increasing N and L, the size of the Floquet matrices scales very unfavor¬ 
ably, i.e. exponentially with these parameters, we must ask ourselves whether it would not be better to use 
a direct propagation method for the given initial state. We use an explicit propagator, e.g. the Runge-Kutta 
(RK) algorithm li45l . More precisely, we resort to a fourth-order RK method with adaptive step-size based 
on a step-to-step error estimation. We applied an error threshold of 10“ 12 for the results shown here, while 
no relevant changes are observed when comparing to thresholds between 1(T 9 • • • 10 -12 . Very importantly, 
we replace the matrix-vector multiplication for the time integration by very efficient vector-vector multipli¬ 
cation. To do this, we store matrices in one-dimensional arrays and just cycle through the upper triangle of 
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Figure 9: (Color online) Polynomial growth of the CPU time with the size of the system d s , for different atom numbers N at fixed 
system sizes L = 5,6 ,..., 12, respectively. For both, the time evolution base on Eq. and RK method the CPU times grow as 
~ d" with exponents n = 3.1 ±0.1 (red line for red/grey symbols) and n = 1.1 ±0.1 (black line for black symbols), respectively. 
The parameters are the same as in Fig. [4] 


the matrix H'(t) in their order of appearance. The corresponding indexes are stored in vectors of integers as¬ 
sociated with the non-zero matrix elements. The number of non-zero elements scales like (~ (8 L + I) x d s ) 
in our case, which saves a lot of allocation memory and considerably reduces the number of operations 
for the time integration. Since the RK algorithm does not intrinsically preserve the norm, we monitor the 
accuracy of the method via the preservation of the norm of the wave function along the time evolution. 

We now compare the efficiency of both propagation methods by evolving in time the upper-band occu¬ 
pation number 


M(t) = (iff 0 \U f (t, 0) ^ ftfUit, 0)|(Ao> = Y J (ils Q \U\t,0)tfU(t,0M Q ) , 


(41) 


l=i 


i=i 


from l = 0 to t ~ 1000 Tn for a fixed Stark force F. The initial condition is given by the lowest energy state 
of the TIF basis for the chosen filling factor N/L, for instance, states of the type |<//o> = |111 ■ ■ ■) ® |000 • • •) 
at filling one. Figure [9] shows that the computation times scale algebraically with d s , as defined around Eq. 
(|9]), for both methods, the Lanczos-based diagonalization with Eq. (37 1 and the RK propagation. The global 
scaling is much better for the RK method for large system sizes, here with exponent n - 1, in contrast to 
the diagonalization with n - 3. As expected, for not too large integration times, the RK method needs less 
CPU time, since here the number of operations is much less than the ones needed for solving the eigenvalue 
problem. Moreover, the storage memory for the RK method is much less than the memory required in the 
Lanczos algorithm, which is already highly optimized with respect to memory requirements. But there is 
a crossover in efficiency measured by the CPU time. The absolute running times scale favorably for the 
Lanczos method only up to system sizes of d s — I (T' 2 , while the RK integration is faster for large sizes 
d s > 10 5 ' 2 , for the fixed overall integration time t = 1000 T « chosen in Fig. [ 9 ] Only for evolutions up to 
very long times, the diagonalization will be favorable again. This is of relevance for small tilts F, since the 
characteristic time scale given by the Bloch period Tg = 1 /F, see also the discussion in section [3] around 
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Figure 10: (Color online) We show good agreement of the two methods detailed in the main text, the spectral (SPA) and long¬ 
time average (LTA) for (a) the upper-band occupation number and (b) the projector operator associated with the state |s ff ) = 
1111 ■ • ■) ® |000 • ■ •). The chosen filling is N/L = 5/5, and we focus on the vicinity of the second order resonance with r = 2. The 
evolution at every value of F is done up to t = 10000 T B . The other parameters are the same as applied in Fig. [4] 


table [T] Of course, the implementation of highly efficient matrix-vector multiplications, as stated above, is 
crucial for the given comparison. We finally note that the observed algebraic scalings of the CPU times are 
stable for a large range of system parameters, which is exemplified as a function of N in Fig. [9] 

Coming back to our discussion in the previous section, we will see that the projection vectors defined 
there are intimately related with the temporal evolution of our system. For our study, we now use the 
Lanczos method which turned out to be more efficient for large evolution times, for which the chaoticity 
is best observed in dynamical signals. The occurrence of avoided crossing is correlated with the strong, 
resonant coupling of at least two quantum states. The evolution of the vector p„ can also be studied by 
means of the projector P a = |.v f> >(xy y |. Based on Eq. (37 1 , we compute the time evolution of the projector 
P„, which reads 


Pa(t) = <«Aol U\t, 0)P a U(t, 0 )|(/r 0 > = |A°(f)| 2 ■ 


(42) 


Flere we set the initial condition to be |</r 0 ) = |s ff ). The system is now let to evolve up to some time /, and 
we compute the long-time average (LTA) 


Pa(t) = 


1 r T 
lim - 

T->oo T J 0 

E W 


PP-hP' 


j+j'k+k 


(43) 


with ojjf = Sj - Ej' and A/^ -k! - k. The integration over the entire time interval basically gives a Dirac 
delta function with argument ojjy - oj^A^k- The delta function will contribute only if ojjf = ANote 
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that A* >k is an integer, which implies that sj = sy + A k'k(*>B maps the eigenenergy sj from the fundamental 
Floquet zone onto a zone at a distance A k'k- The Floquet zones are equivalent, then without loss of generality 
Sj — » sy. Since degenerancies arc excluded in our spectra due to just ACs in the chaotic regime, the second 
right-hand term in Eq. ( [43] ) vanishes in the LTA, and we obtain 

\cf\dj -> 2 p a j(sj\P a \sj) = ipr a (F) . (44) 

jk j 


This means that the system becomes ergodic during the evolution [46 ]. The effective number of participating 
states is hence large, and, in the extreme case, the vector pj = (• • • ,p°j, ■ ■ -j approaches the equipartition 


condition p' = d~ l (1,1, • • •, 1, • • •, 1,1). We conclude that both, the appearance of ACs in the spectrum 
and the spectral averages (SPA) of P a = |.v (( ,)(y t .| reflect themselves directly in the long-time behavior of the 
projection vectors. 

The LTA analysis can be straightforwardly extended to any other quantity, for instance the occupation 
number M of the upper energy band, whose long-time average reduces to M(t) ~ J/j P a (£j\M\sj). Figure 


10 shows the comparison of the actual values and the ergodic limits of the observables mentioned above 


for N/L = 5/5 (d s = 402), around a second order resonance for initial state lAo) = 1111 • • •) <8> |000 - • ■). 
In order to numerically compute M{t) and P a (t) we evolve the initial state up to time t = 10000Tb. Note 
that for both, M(t ) and P„(t) their LTA and SPA correspond best close to ACs, that is whenever 1 /ipr ff is 
maximal. Linally, by evolving a given initial state \i//q) it is possible to detect the signatures of avoided 
crossings within the interesting spectral regions (here close to RET conditions). This gives an alternative 
and most likely easier route to the experimental measurements of signature of ACs and hence of quantum 
chaos in complex quantum systems, in the same spirit as the behavior of Bloch oscillation does in simpler 
single-band systems fR .iTBl. 


6. Conclusions 

In this paper we have presented exact numerical methods for the treatment of a one-dimensional 
Wannier-Stark system including strong interactions and two Bloch bands. This system has been of great ex¬ 
perimental interest over the last decade iTTl l9l ITTl IT2l IT3l l47ll . and experiments studying the strong coupling 
in the chaotic regime arc just on their way llT4l . Lor this implicitly time-dependent quantum system (since 
driven system by the Stark force F), we implemented two independent numerical approaches, one based 
on the combination of the Lloquet theory and the Lanczos algorithm for symmetric matrices, and a second 
one based on a Runge-Kutta integration of a given initial state. In both cases, we ensured efficient matrix- 
vector multiplications and the memory storage of only non-zero elements in the corresponding Hamiltonian 
matrices. The computational times scale polynomially with the system size in both cases. Depending on 
the dimension of the system one or the other algorithm is more efficient for the temporal evolution over 
reasonable integration times. Lor large system sizes and not extremely long integration times, the direct and 
explicit Runge-Kutta integration proved to be more efficient, scaling with an exponent of approximately 
one with the system size. Nevertheless, the spectral analysis of the Lloquet Hamiltonian is very useful for 
characterizing the dynamical properties. Interestingly, signatures of the chaoticity of the system can also be 
found in the temporal evolution of appropriate observables. The avoided crossing scenarios in the spectrum 
as a function of the control parameter F arc directly related to the dynamical spreading across a given basis. 
This spreading even leads to complete ergodicity in the quantum chaotic regime of our model occurring at 
resonant tunneling conditions between the two energy bands. 

The numerical methods presented here are suitable for a larger variety of problems similar to the multi¬ 
band version of the Wannier-Stark system. Lor instance, periodic driving of bosonic many-body systems 
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could be used to create highly entangled states f48l . Moreover, artificial gauge fields can be created by Stark 
forces or by time-periodic potentials |[49ll in one or more spatial dimensions, possibly also in the strongly 
interacting regime with comparable resonant forces in the near future.Generally, as we showed for our 
case, the numerical performance may be improved exploiting the characteristic energy scales of the given 
system for an optimal choice of the numerical parameters. Our work shows that generic quantum lattice 
systems with strong many-body interactions should be tractable by exact methods up to Hilbert space sites 
of about 10 5 on standard scientific workstations. Using sophisticated parallelization techniques and larger 
computing clusters with access to much more memory, one might gain at least another order of magnitude 
for the maximally treatable Hilbert space sizes Il30l . 
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